
library(ggplot2)
require(data.table)
library("dplyr")
library("ggpubr") 
require(tidyverse) #for easy data manipulation and visualization
require(caret) #for easy machine learning workflow
#library(ncdf4)

wd <- "C:\\Users\\YFan\\...\\NFood_replication\\"

pft = c(17,18,  19,20,   23,24,    41,42,    61,62,   67,68,   75,76, 77,78)#rainfed vs. irrigated
crops = c("Corn", "Wheat", "Soybean", "Cotton", "Rice", "Sugarcane", "Tropical Corn", "Tropical Soybean")
cropname = c("Maize", "Wheat", "Soy", "Cotton", "Rice", "Sugarcane")

case = c("RCP45","RCP45_85LU","RCP85","RCP85_45CO2","RCP85_SAI","RCP85_MSB","RCP85_CCT")
casename = c("RCP45",expression('RCP45'[SSP5]),"RCP85",expression('RCP85'[RCP45CO2]),"RCP85+SAI","RCP85+MSB","RCP85+CCT")


#******************************************************************************
#=============================Statistics and Figures===========================
#******************************************************************************

load(paste0(wd,"/data/RData/Fig2.Crop_yield_validation_FAO.RData"))

df0= rbind(md.yield,fao.yield) #for facet plot
df1= fao.yield[md.yield, on=.(Crop,Year)] #join two tables, for scatterplot


#==================plot option 1==================

#####Production and Area validation (draw individual crops one by one with consistent scales on production and area)
#===identical scale for all crops
svg(paste0(wd,"/figures/main/Fig.2a.Production_Area_FAO_validation.svg"), width = 5.8, height = 5.)
#jpeg(paste0(wd,"/figures/revision/Fig.2a.Production_Area_FAO_validation.jpeg"), units = "in", width =6, height=10, res = 600)

split.screen(c(3,2))
screen(1) # prepare screen 1 for output
par(mar=c(.5, 3.5, 2, 1.8) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9, col="black")
plot(TotProd~Year, data=df0[Crop=="Maize"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Maize", cex.main=1, font.main=1,
     ylim=c(100, 1100), xlab="", ylab="") 
title("Maize", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(600,1000,200))
lines(TotProd~Year, data=df0[Crop=="Maize"& Case=="CLM5"], lty="dashed") #
#legend("center", c("FAO", "CLM5"), lty=c("solid", "dashed"), bty = "n", horiz = TRUE)
par(new = TRUE, xpd = TRUE, col="#D55E00")
plot(TotArea~Year, data=df0[Crop=="Maize"& Case=="FAO"], type="l", bty="n", xaxt="n", yaxt="n", 
     ylim=c(100,400), xlab="", ylab="")
lines(TotArea~Year, data=df0[Crop=="Maize"& Case=="CLM5"], lty="dashed") #
axis(side = 4, at=seq(100,200,50), col.ticks ="#D55E00", col.axis="#D55E00")

screen(2) # prepare screen 1 for output
par(mar=c(.5, 1.8, 2, 3.5) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9, col="black")
plot(TotProd~Year, data=df0[Crop=="Sugarcane"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Sugarcane", cex.main=1, font.main=1,
     ylim=c(-500, 500), xlab="", ylab="") 
title("Sugarcane", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(0,400,200))
lines(TotProd~Year, data=df0[Crop=="Sugarcane"& Case=="CLM5"], lty="dashed") #
legend("top", c("FAO", "CLM5"), lty=c("solid", "dashed"), bty = "n", cex=0.9, horiz = TRUE)
par(new = TRUE, xpd = TRUE, col="#D55E00")
plot(TotArea~Year, data=df0[Crop=="Sugarcane"& Case=="FAO"], type="l", bty="n", xaxt="n", yaxt="n", 
     ylim=c(0,300), xlab="", ylab="")
lines(TotArea~Year, data=df0[Crop=="Sugarcane"& Case=="CLM5"], lty="dashed") #
axis(side = 4, at=seq(0,100,50), col.ticks ="#D55E00", col.axis="#D55E00")

screen(3) # prepare screen 1 for output
par(mar=c(1., 3.5, 1.5, 1.8) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9, col="black", xpd = TRUE)
plot(TotProd~Year, data=df0[Crop=="Wheat"& Case=="FAO"], type="l", xaxt="n", yaxt="n",# main="Wheat", cex.main=1, font.main=1,
     ylim=c(-100, 900), xlab="", ylab="")
title(main="Wheat", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(400,800,200))
lines(TotProd~Year, data=df0[Crop=="Wheat"& Case=="CLM5"], lty="dashed") #
par(new = TRUE, xpd = TRUE, col="#D55E00")
plot(TotArea~Year, data=df0[Crop=="Wheat"& Case=="FAO"], type="l", bty="n", xaxt="n", yaxt="n", 
     ylim=c(200,500), xlab="", ylab="")
lines(TotArea~Year, data=df0[Crop=="Wheat"& Case=="CLM5"], lty="dashed") #
axis(side = 4, at=seq(200,300,50), col.ticks ="#D55E00", col.axis="#D55E00")

screen(4) # prepare screen 1 for output
par(mar=c(1., 1.8, 1.5, 3.5) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9, col="black")
plot(TotProd~Year, data=df0[Crop=="Rice"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Rice", cex.main=1, font.main=1,
     ylim=c(100,1100), xlab="", ylab="") 
title(main="Rice", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(600,1000,200))
lines(TotProd~Year, data=df0[Crop=="Rice"& Case=="CLM5"], lty="dashed") #
par(new = TRUE, xpd = TRUE, col="#D55E00")
plot(TotArea~Year, data=df0[Crop=="Rice"& Case=="FAO"], type="l", bty="n", xaxt="n", yaxt="n", 
     ylim=c(100,400), xlab="", ylab="")
lines(TotArea~Year, data=df0[Crop=="Rice"& Case=="CLM5"], lty="dashed") #
axis(side = 4, at=seq(100,200,50), col.ticks ="#D55E00", col.axis="#D55E00")
corners = par("usr") #Gets the four corners of plot area (x1, x2, y1, y2)
text(x = corners[2]+3.5, y = mean(corners[3:4]), "Area (million ha)", col ="#D55E00", srt = 270) #mtext does not accept "srt"

screen(5) # prepare screen 1 for output
par(mar=c(1.5, 3.5, 1, 1.8) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9, col="black")
plot(TotProd~Year, data=df0[Crop=="Soy"& Case=="FAO"], type="l", xaxt="n", yaxt="n", # main="Soy", cex.main=1, font.main=1,
     ylim=c(-500, 500), xlab="", ylab="") 
title(main="Soy", line=0.5, cex.main=1, font.main=1)
axis(side = 1, at=c(2006,2012,2018))
axis(side = 2, at=seq(0,400,200))
lines(TotProd~Year, data=df0[Crop=="Soy"& Case=="CLM5"], lty="dashed") #
par(new = TRUE, xpd = TRUE, col="#D55E00")
plot(TotArea~Year, data=df0[Crop=="Soy"& Case=="FAO"], type="l", bty="n", xaxt="n", yaxt="n", 
     ylim=c(50,350), xlab="", ylab="")
lines(TotArea~Year, data=df0[Crop=="Soy"& Case=="CLM5"], lty="dashed") #
axis(side = 4, at=seq(50,150,50), col.ticks ="#D55E00", col.axis="#D55E00")

screen(6) # prepare screen 1 for output
par(mar=c(1.5, 1.8, 1, 3.5) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9, col="black")
plot(TotProd~Year, data=df0[Crop=="Cotton"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Cotton", cex.main=1, font.main=1,
     ylim=c(-500, 500), xlab="", ylab="")
title(main="Cotton", line=0.5, cex.main=1, font.main=1)
axis(side = 1, at=c(2006,2012,2018))
axis(side = 2, at=seq(0,400,200))
lines(TotProd~Year, data=df0[Crop=="Cotton"& Case=="CLM5"], lty="dashed") #
# legend("top", c("FAO", "CLM5"), lty=c("solid", "dashed"), bty = "n", horiz = FALSE)
par(new = TRUE, xpd = TRUE, col="#D55E00")
plot(TotArea~Year, data=df0[Crop=="Cotton"& Case=="FAO"], type="l", bty="n", xaxt="n", yaxt="n", 
     ylim=c(0,300), xlab="", ylab="")
lines(TotArea~Year, data=df0[Crop=="Cotton"& Case=="CLM5"], lty="dashed") #
axis(side = 4, at=seq(0,100,50), col.ticks ="#D55E00", col.axis="#D55E00")

close.screen(all = TRUE)

par(mar=c(0, 2, 1.5, 1), mgp=c(0.5,0.5,0), las=1, new=TRUE) #overlay plot to the above split plot
plot(TotProd~Year, data=df0, type="n", bty="n", xaxt="n", yaxt="n", xlab="", ylab=expression("Production (million tonnes year "^-1*")"), 
     main=expression('('*bold(a)*')       '), cex=0.9)

dev.off()



#===identical scale for all crops
svg(paste0(wd,"/figures/main/Fig.2b.AvgYield_FAO_validation.svg"), width = 5., height = 5.)

split.screen(c(3,2))
screen(1) # prepare screen 1 for output
par(mar=c(.5, 3, 2, 0) + 0.1, mgp=c(1.5,0.5,0), las=1, cex =0.9, col="black")
plot(AvgYield~Year, data=df0[Crop=="Maize"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Maize", cex.main=1, font.main=1,
     ylim=c(3.8, 5.8), xlab="", ylab="") 
title("Maize", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(4,5.5,0.5))
lines(AvgYield~Year, data=df0[Crop=="Maize"& Case=="CLM5"], lty="dashed") #
legend("top", c("FAO", "CLM5"), lty=c("solid", "dashed"), bty = "n", cex=0.9, horiz = TRUE)

screen(2) # prepare screen 1 for output
par(mar=c(.5, 2, 2, 1) + 0.1, mgp=c(1.5,0.5,0), las=1, cex =0.9, col="black")
plot(AvgYield~Year, data=df0[Crop=="Sugarcane"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Sugarcane", cex.main=1, font.main=1,
     ylim=c(5.2, 7.2), xlab="", ylab="") 
title("Sugarcane", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(5.5,7,0.5))
lines(AvgYield~Year, data=df0[Crop=="Sugarcane"& Case=="CLM5"], lty="dashed") #

screen(3) # prepare screen 1 for output
par(mar=c(1., 3, 1.5, 0) + 0.1, mgp=c(1.5,0.5,0), las=1, cex =0.9, col="black", xpd = TRUE)
plot(AvgYield~Year, data=df0[Crop=="Wheat"& Case=="FAO"], type="l", xaxt="n", yaxt="n",# main="Wheat", cex.main=1, font.main=1,
     ylim=c(1.8, 3.8), xlab="", ylab="")
title(main="Wheat", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(2,3.5,0.5))
lines(AvgYield~Year, data=df0[Crop=="Wheat"& Case=="CLM5"], lty="dashed") #

screen(4) # prepare screen 1 for output
par(mar=c(1., 2, 1.5, 1) + 0.1, mgp=c(1.5,0.5,0), las=1, cex =0.9, col="black")
plot(AvgYield~Year, data=df0[Crop=="Rice"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Rice", cex.main=1, font.main=1,
     ylim=c(3.7, 5.7), xlab="", ylab="") 
title(main="Rice", line=0.5, cex.main=1, font.main=1)
axis(side = 2, at=seq(4,5.5,0.5))
lines(AvgYield~Year, data=df0[Crop=="Rice"& Case=="CLM5"], lty="dashed") #

screen(5) # prepare screen 1 for output
par(mar=c(1.5, 3, 1, 0) + 0.1, mgp=c(1.5,0.5,0), las=1, cex =0.9, col="black")
plot(AvgYield~Year, data=df0[Crop=="Soy"& Case=="FAO"], type="l", xaxt="n", yaxt="n", # main="Soy", cex.main=1, font.main=1,
     ylim=c(1.8, 3.8), xlab="", ylab="") 
title(main="Soy", line=0.5, cex.main=1, font.main=1)
axis(side = 1, at=c(2006,2012,2018))
axis(side = 2, at=seq(2,3.5,0.5))
lines(AvgYield~Year, data=df0[Crop=="Soy"& Case=="CLM5"], lty="dashed") #

screen(6) # prepare screen 1 for output
par(mar=c(1.5, 2, 1, 1) + 0.1, mgp=c(1.5,0.5,0), las=1, cex =0.9, col="black")
plot(AvgYield~Year, data=df0[Crop=="Cotton"& Case=="FAO"], type="l", xaxt="n", yaxt="n", #main="Cotton", cex.main=1, font.main=1,
     ylim=c(-0.2, 1.8), xlab="", ylab="")
title(main="Cotton", line=0.5, cex.main=1, font.main=1)
axis(side = 1, at=c(2006,2012,2018))
axis(side = 2, at=seq(0.0,1.5,0.5))
lines(AvgYield~Year, data=df0[Crop=="Cotton"& Case=="CLM5"], lty="dashed") #

close.screen(all = TRUE)

par(mar=c(0, 2, 1.5, 0.5), mgp=c(0.5,0.5,0), las=1, new=TRUE) #overlay plot to the above split plot
plot(AvgYield~Year, data=df0, type="n", bty="n", xaxt="n", yaxt="n", xlab="", ylab=expression("Yield (tonnes ha"^-1*")"), 
     main=expression('('*bold(b)*')'), cex=0.9)
#legend("bottom", c("FAO", "CLM5"), lty=c("solid", "dashed"), bty = "n", horiz = TRUE)
dev.off()





# #==================plot option 2==================
# #========ggplot: hard to make uniform y scales 
# #In order to change the order of the panels, you need to change the factor levels for your faceting variable Crop
# cropname <- c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")
# 
# left_breaks <- function(x) { if (max(x) > 500) seq(500, 1000, 250) else if (max(x) > 200) seq(200, 300, 50) else if (max(x) > 100) seq(100, 160, 60) else seq(15, 25, 10) }
# right_breaks <- function(x) {if (max(x) > 600) seq(150, 250, 100) else if (max(x) > 400) seq(200, 300, 100) else if (max(x) > 200) seq(90, 120, 30)  else seq(20, 35, 15) }
# #right breaks also according to level of TotProd
# df0[Year==2018,.(Case,Crop,TotProd)] 
# 
# #df<- df0[, TotArea:= TotArea*2] #scaled for sec_axis
# #use mutate to internally reorder Crop based on specified levs
# p0 <-  df0 %>% 
#   mutate(Crop = factor(Crop, levels=cropname, labels=cropname)) %>%
#   mutate(Case = factor(Case, levels=c("FAO", "CLM5"))) %>%
#   ggplot(aes(Year, TotProd, group=Case)) + geom_line(aes(linetype = Case), color = "black", size = 0.7)+ 
#   #change linetype, color, size mannually
#   scale_linetype_manual(name = "Case",  values = c("solid", "dashed"))  +
#   #additonal scale for crop area
#   geom_line(aes(Year, TotArea*1.6, group=Case, linetype=Case), color = "#D55E00", size = 0.7 ) + 
#   facet_wrap(~Crop, nrow=3, scales='free_y') +  
#   #facet_grid_sc(rows=vars(Crop), scales = list(y = scales_y_left)) + #set different scales for subplots
#   scale_x_continuous(breaks =c(2006,2012,2018)) +
#   scale_y_continuous(expression("Production (million tonnes year "^-1*")"), breaks=left_breaks,
#                      expand = expand_scale(mult = c(0.3, 0.1), add = c(0, 0)),
#                      sec.axis = sec_axis(~ . *1/1.6, name = "Area (million ha)",breaks=right_breaks))
# #merge the legends for color and linetype 
# p0 <- p0 + guides(linetype = guide_legend( keywidth = 2, label.position = "right", label.hjust = 0, 
#                                            override.aes = list(size = 0.7, color="black", linetype = c("solid", "dashed")) )) 
# ##remove grid and background or remove all
# p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) 
# #p <- p + xlab("Year") + ylab("Yield change (MtC/year)")
# p0 <- p0 + ggtitle(expression('('*bold(a)*')')) +
#   theme(strip.background = element_blank(), strip.text = element_text(size = 12))+
#   theme(plot.title = element_text(size = 16, hjust=0.5), panel.spacing = unit(0, "lines")) +
#   theme(legend.text = element_text(size=15), legend.title = element_blank(), legend.position = "bottom") + 
#   theme(axis.text = element_text(size=12), axis.title=element_text(size=14))
# # Change the color palette
# p0 +   theme(axis.title.y.left = element_text(color = "black"),
#              axis.text.y.left = element_text(color = "black"),
#              axis.title.y.right = element_text(color = "#D55E00"), axis.ticks.y.right = element_line(color = "#D55E00"),
#              axis.text.y.right = element_text(color = "#D55E00"))
# 
# #or use ggsave to save the last plot on the screen
# ggsave(paste0(wd,"/crop_figures/revision/Fig.2a.Production_Area_FAO_validation.svg"), device ="svg", units = "in", width = 7, height = 6, dpi = 600)
# #ggsave(paste0(getwd(),"/Manucript1/Figures/Fig.Ext1a.Production_Area_FAO_validation.eps"), device ="eps", units = "in", width = 7, height = 6, dpi = 600)
# 
# 
# 
# #=====linegraph of global crop yields modeled vs. FAO data
# #jpeg(paste0(getwd(),"/crop_figures/final/Fig.Ext1.AvgYield_FAO_validation.jpeg"), units = "in", width = 6, height = 6, res = 600)
# 
# #In order to change the order of the panels, you need to change the factor levels for your faceting variable Crop
# cropname <- c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")
# 
# #use mutate to internally reorder Crop based on specified levs
# p0 <-  df0 %>% 
#   mutate(Crop = factor(Crop, levels=cropname, labels=cropname)) %>%
#   mutate(Case = factor(Case, levels=c("FAO", "CLM5") )) %>%
#   ggplot(aes(Year, AvgYield, group=Case)) + geom_line(aes(linetype = Case), color = "black", size = 0.7)+ 
#   #change linetype, color, size mannually
#   scale_linetype_manual(name = "Case",  values = c("solid", "dashed"))  +
#   facet_wrap(~Crop, nrow=3, scales='free_y') +  
#   #facet_wrap(~Crop, nrow=3, ncol=2, scales='fixed') + 
#   scale_x_continuous(name="",breaks =c(2006,2012,2018)) +
#   scale_y_continuous(expression("Yield (tonnes ha"^-1*")"), labels = function(x) sprintf("%.1f", x))
# #merge the legends for color and size 
# p0 <- p0 + guides(linetype = guide_legend( keywidth = 2, label.position = "right", label.hjust = 0, 
#                                            override.aes = list(size = 0.7, color="black", linetype = c("solid", "dashed")) )) 
# #override linetype in previous guides for yield and area
# ##remove grid and background or remove all
# p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) 
# #p <- p + xlab("Year") + ylab("Yield change (MtC/year)")
# p0 <- p0 + ggtitle(expression('('*bold(b)*')')) +
#   theme(strip.background = element_blank(), strip.text = element_text(size = 12))+
#   theme(plot.title = element_text(size = 16,hjust=0.5), panel.spacing = unit(1, "lines")) +
#   theme(axis.text = element_text(size=12), axis.title=element_text(size=14)) +
#   #theme(legend.text = element_text(size=15), legend.title = element_blank(), legend.position = "bottom") 
#   theme(legend.text = element_text(size=15), legend.title = element_blank(), legend.position = c(0.8, 0.15), legend.background = element_blank()) 
# 
# # Change the color palette
# #p0 <- set_palette(p0, palette = cbPalette[-1])
# p0 +   theme(axis.title.y.left = element_text(color = "black"),
#              axis.text.y.left = element_text(color = "black"))
# 
# #dev.off()
# ggsave(paste0(wd,"/crop_figures/revision/Fig.2b.AvgYield_FAO_validation.svg"), device ="svg", units = "in", width = 6, height = 6, dpi = 600)
# #ggsave(paste0(getwd(),"/Manucript1/Figures/Fig.Ext1b.AvgYield_FAO_validation.eps"), device ="eps", units = "in", width = 6, height = 6, dpi = 600)
# 
# 
# 
# #==================plot option 3==================
# #=======draw all crops in one plot with different colors or linetypes to save space
# 
# library(viridis)
# #======final color for crops=====
# cropcolor <- viridis(18)[seq(7,18,2)]
# barplot(height=rep(1,length(cropcolor)),col=cropcolor) 
# 
# cropname <- c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")
# 
# #==Production==
# svg(paste0(wd,"/crop_figures/revision/Fig.2a.Production_FAO_validation-v3.svg"), width = 3.5, height = 5)
# 
# p0 <-  df0 %>% 
#   mutate(Crop = factor(Crop, levels=cropname, labels=cropname)) %>%
#   mutate(Case = factor(Case, levels=c("FAO", "CLM5"))) %>%
#   ggplot(aes(Year, TotProd)) + geom_line(aes(linetype = Case, color = Crop), size = 0.7)+ geom_point(aes(shape = Crop, color = Crop))+
#   #change linetype, color, size mannually
#   scale_color_manual(values = cropcolor)  +
#   scale_shape_manual(values=c(0,8,1,17,16,2)) +
#   scale_linetype_manual(name = "Case",  values = c("solid", "dashed"))  +
#   # #additonal scale for crop area
#   # geom_line(aes(Year, (TotArea-400)*2.5, color = Crop, linetype=Case), alpha = 0.3, size = 0.8 ) + 
#   # geom_point(aes(Year, (TotArea-400)*2.5, shape = Crop, color = Crop)) +
#   scale_x_continuous(name="", breaks =c(2006,2012,2018)) +
#   scale_y_continuous(expression("Production (million tonnes year "^-1*")"), limit=c(0,1100), breaks=seq(0,1000,200))  #left axis breaks
#                      #expand = expand_scale(mult = c(0.3, 0.1), add = c(0, 0)),
#                      #sec.axis = sec_axis(~ ./2.5+400, name = "Area (million ha)", breaks=seq(0,300,50)))
# #merge the legends for color and linetype 
# p0 <- p0 + guides(linetype = guide_legend( keywidth = 2, label.position = "right", label.hjust = 0, ncol=2),
#                   color = guide_legend( keywidth = 2, keyheight=0.7, label.position = "right", label.hjust = 0, ncol=2))
# 
# ##remove grid and background or remove all
# p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) 
# 
# p0 + ggtitle(expression('('*bold(a)*')')) +
#   theme(strip.background = element_blank(), strip.text = element_text(size = 12))+
#   theme(plot.title = element_text(size = 14, hjust=0.5), panel.spacing = unit(0, "lines")) +
#   theme(axis.text = element_text(size=12), axis.title=element_text(size=14)) +
#   theme(legend.text = element_text(size=10), legend.title = element_blank(), legend.background = element_blank(),
#         legend.position = c(0.45, 0.78), legend.spacing.y = unit(2.2, "cm"))  #adjust the spacing of two different legends
#   
# dev.off()
# 
# 
# #==Area==
# svg(paste0(wd,"/crop_figures/revision/Fig.2b.CropArea_FAO_validation-v3.svg"), width = 3.5, height = 5)
# p0 <-  df0 %>% 
#   mutate(Crop = factor(Crop, levels=cropname, labels=cropname)) %>%
#   mutate(Case = factor(Case, levels=c("FAO", "CLM5") )) %>%
#   ggplot(aes(Year, TotArea)) + geom_line(aes(linetype = Case, color = Crop), size = 0.8)+ geom_point(aes(shape = Crop, color = Crop))+
#   #change linetype, color, size mannually
#   scale_color_manual(values = cropcolor)  +
#   scale_shape_manual(values=c(0,8,1,17,16,2)) +
#   scale_linetype_manual(name = "Case",  values = c("solid", "dashed"))  +
#   scale_x_continuous(name="",breaks =c(2006,2012,2018)) +
#   scale_y_continuous(expression("Area (million ha)"), limits=c(0,300), breaks=seq(0,300,50), labels = function(x) sprintf("%.0f", x))
# #merge the legends for color and size 
# p0 <- p0 + guides(linetype = guide_legend( keywidth = 3, label.position = "right", label.hjust = 0))
# 
# ##remove grid and background or remove all
# p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) 
# #p <- p + xlab("Year") + ylab("Yield change (MtC/year)")
# p0 + ggtitle(expression('('*bold(b)*')')) +
#   theme(strip.background = element_blank(), strip.text = element_text(size = 12))+
#   theme(plot.title = element_text(size = 14,hjust=0.5), panel.spacing = unit(1, "lines")) +
#   theme(axis.text = element_text(size=12), axis.title=element_text(size=14)) +
#   theme(legend.text = element_text(size=10), legend.title = element_blank(), legend.position = "none") 
# #theme(legend.text = element_text(size=10), legend.title = element_blank(), legend.position = c(0.8, 0.15), legend.background = element_blank()) 
# dev.off()
# 
# 
# 
# 
# #==Yield==
# svg(paste0(wd,"/crop_figures/revision/Fig.2c.AvgYield_FAO_validation-v3.svg"), width = 3.5, height = 5)
# cropname <- c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")
# #use mutate to internally reorder Crop based on specified levs
# p0 <-  df0 %>% 
#   mutate(Crop = factor(Crop, levels=cropname, labels=cropname)) %>%
#   mutate(Case = factor(Case, levels=c("FAO", "CLM5") )) %>%
#   ggplot(aes(Year, AvgYield)) + geom_line(aes(linetype = Case, color = Crop), size = 0.8)+ geom_point(aes(shape = Crop, color = Crop))+
#   #change linetype, color, size mannually
#   scale_color_manual(values = cropcolor)  +
#   scale_shape_manual(values=c(0,8,1,17,16,2)) +
#   scale_linetype_manual(name = "Case",  values = c("solid", "dashed"))  +
#   scale_x_continuous(name="",breaks =c(2006,2012,2018)) +
#   scale_y_continuous(expression("Yield (tonnes ha"^-1*")"), labels = function(x) sprintf("%.1f", x))
# #merge the legends for color and size 
# p0 <- p0 + guides(linetype = guide_legend( keywidth = 3, label.position = "right", label.hjust = 0, 
#                                            override.aes = list(size = 0.8, color="black", linetype = c("solid", "dashed")) )) 
# #override linetype in previous guides for yield and area
# ##remove grid and background or remove all
# p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) 
# #p <- p + xlab("Year") + ylab("Yield change (MtC/year)")
# p0 + ggtitle(expression('('*bold(c)*')')) +
#   theme(strip.background = element_blank(), strip.text = element_text(size = 12))+
#   theme(plot.title = element_text(size = 14,hjust=0.5), panel.spacing = unit(1, "lines")) +
#   theme(axis.text = element_text(size=12), axis.title=element_text(size=14)) +
#   theme(legend.text = element_text(size=10), legend.title = element_blank(), legend.position = "none") 
#   #theme(legend.text = element_text(size=15), legend.title = element_blank(), legend.position = c(0.8, 0.15), legend.background = element_blank()) 
# dev.off()



